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In this paper we theoretically study propagation of steady state ultrashort pulse in dissipative 
medium. We considered two cases (i) medium consists of lossy metallic nanostructures embedded 
into a gain material and (ii) the gain material is embedded directly into the nanostructures. We 
found the shape and velocity of an optical pulse coupled with the polarization wave. 

PACS numbers: 



INTRODUCTION 

For over a decade extensive research has been undertaken 
on negative refractive indexed metamaterials [1, 2, 3] 
both in the microwave and optical regimes. The research 
in optical metamaterials has been stimulated by the pos- 
sibility of light manipulation methods and potential ap- 
plications. Most of the experimental work with metama- 
terials utilizes the plasmonic resonance in metallic struc- 
tures embedded into a host dielectric. The disadvantage 
of this type of metamaterial arises from the high loss 
value that results from an induced plasmonic resonance 
in the optical range. Currently, the reduction or com- 
pensation of these losses has been the the focus of inten- 
sive research. Conceptually, the compensation of losses 
can be achieved in two ways. First, by adding gain into 
a host material containing the metallic nanostructures. 
Second, by adding a gain material into the nanostruc- 
tures themselves. In this paper, we consider the optical 
pulse dynamics in a negative index metamaterial using 
both types of loss compensation techniques. To under- 
stand the general features of the pulse dynamics in such 
a system we consider a simplified mathematical model. 



BASIC EQUATIONS 

The propagation of an electromagnetic wave in matter 
is described by the system of Maxwell's equations: 
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1 dB 
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By choosing a system of coordinates in such a way as 
to assume that a wave propagates in the z-direction, 
the magnetic field is oriented in y-direction, H = 
(0, H(z, t), 0), and the electric field is oriented in the re- 
direction, E — (E(z, t), 0, 0), we transform the system of 
vector equations (1) into a system of scalar equations: 
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(2) 



We start our analysis of the system from the initial case 
of an electromagnetic wave as it propagates through a 
nanostructured host dielectric without a gain material. 
We also assume that all resonances of the host material 
are far from the resonances of the metallic nanostruc- 
tures. In this case, D and B can be represented in the 
standard form through the medium polarization, P, and 
the magnetization, M: 



B = H + AttM , D = E + AnP. 



(3) 



To close the system of equations we need to add ma- 
terial equations that describe the interaction of P and 
M with an electromagnetic wave. The total medium re- 
sponse can be represented as sum of the responses from 
the nanostructures and the gain material: P = P ns + P g 
and M = M ns . To describe such an interaction, we 
consider the simplified model represented by distributed 
linear electric and magnetic resonators. The plasmonic 
oscillations result in an oscillation of the polarization de- 
scribed by following equation [4]: 
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Here to p is the effective plasma frequency, u>d is the fre- 
quency of the dimensional quantization that results from 
the geometry of the resonator, and 7 accounts for losses. 
In the popular fishnet structures the dielectric permittiv- 
ity is negative in operating frequency range. In this case 
we can choose u>d = 0. 

To describe the magnetic response we assume that a 
magnetic resonator in a leading order operates as an LC- 
circuit [5]. This allows us to write a partial differential 
equation that relates the magnetization and the magnetic 
field as follows [61: 



can therefore be expressed 



d 2 M n 

dt 2 



0M n 

( -dt 



+ u? T M ns = -b 



d 2 H 



(5) 



where L, C are effective inductance and capacity, lo\ is a 
Tompson frequency, a is a dissipative constant, and b de- 
scribes coupling of nanoresonator with external magnetic 
flux [7]. 

The system (2), (3), (4), (5) is a generalization of the 
well known Maxwell-Lorentz system of equations, which 
describes the resonant interaction of an electromagnetic 
field with a system of linear oscillators that responds only 
to the electric component of incident light. This system 
has played a fundamental role in understanding the in- 
teraction of light with resonant atoms [10, 11]. 

The host material is characterized by the dielectric per- 
mittivity, £o(w), and the magnetic permeability, fio(u)). 
For most materials the permeability in the optical range 
is 1, that is, fio(w) = 1. Now, by taking the Fourier 
transform of the equations in (2) and substituting in the 
expressions for B and D we obtain the following equa- 
tions for Fourier components E and H in terms of P and 
M: 



k 2 E 



k 2 H 
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-kP. (6) 



In order to describe the response from the nanostruc- 
tures we introduce both an effective dielectric and a mag- 
netic susceptibility. We consider here the case of a "di- 
luted" metamaterial when near neighbor interactions of 
the nanostructurcs arc much weaker than their interac- 
tion via the electromagnetic field. The nanostructure 
contribution to the effective permittivity and permeabil- 
ity can be obtained from equations (4) and (5) by apply- 
ing the Fourier transform, that is, 



P — 



(ujjj — LU 2 ) — lLJ"f 
M ^ 



E = Xe (uj)E, 



H = Xm {u)H. 



(7) 



(8) 



e{uj) = £ +47TXeM, 
jl(u>) = l+47TXm(w) 



(9) 



Using (9), the system of equations (6) can be represented 
as 



k 2 



E(k,i 



47T 
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H(k,(j) = 



4:TT — kPg(k,LU). 



(11) 



Ultrashort optical pulses are well described by a slowly 
varying envelope approximation. This approximation al- 
lows one to simplify the system (2), (3), (4), (5) by 
representing E, H , P and M in the following form: 

E(z,t) = £(z 1 1) exp (iujot — ik^z) + ex., 

P(z,t) = V(z, t) exp (iuot — ikoz) + c.c, 

H(z,t) = H(z,t)exp(iuJot — ikoz) + c.c., (12) 

M(z,t) — M(z,t)exp(iuot — ikoz) + c.c., 

where £(z,t), V(z,t), H(z,t) and M(z,t) are slowly- 
varying functions of t and z. It should be noted that the 
constant ko stands for the projection of the wave vector 
on z axis, and that it can be both positive and negative 
for positive and negative index materials respectively. 

The factorization (12) means that the spectral width 
of the signal is much smaller than carrier frequency, 
therefore the equations for the slowly varying amplitudes 
can be obtained by a frequency and wave vector shift, 
ui — > u + ujq, k — > k + k, in Fourier space [12]. Thus, 
the equations describing the slow dynamics read 



k 2 



£(k, ui) 



4^(^) 2 jHu)V g {k,u) (13) 



- ) e{ib)jl{Cb) 



H(k,u>) 



The effective permittivity, e, and the permeability, ft, 



A-K^kV g {k,uj), (14) 

Cj = u> + loq, k = k + ko- 

Here the parameter fco can be both positive and negative, 
since it represents the projection of the wave vector of the 
carry wave on the direction of the wave propagation. In 
the case of a negative refractive index the sign of ko is 
negative. 

Since the width of wave packet is much less than carrier 
frequency, the inequalities lo <C loq and k -c |fco| still hold. 
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By expanding the expression in the square brackets in 
equations (13) we obtain 



/cq -\- 2kok T k 2 



(u>l + 2lu q lo + J 2 ) 



n 2 (uj ) 



+ 



dn 2 
div 



1 



du 2 



(15) 



where n 2 = e(u))jl(uj) is a square of the effective complex 
index of refraction. 

The effective index of refraction has both a real and an 
imaginary part. The imaginary part is responsible for the 
phase dynamics of the field, while the imaginary part is 
responsible for field decay. We consider the propagation 
constant fc, which is defined as follows 



(16) 



k(ui) — (u}/c)Re^/e(co)ij(ui). 

Introducing the auxiliary function T(uj) 

T(uj) — (w/c)Im\/£(w)p(w) 

and applying the inverse Fourier transform, we obtain 
equations in spatio-temporal variables: 




d 2 k 



d 2 
dt 2 



27rwoA i (^o) 




c Re y/(i(u} Q )e{u ) 

d 2 k 
duo 2 



V, 



9> 



(17) 




Here the group velocity v g is defined as v^ 1 - 
and To = T(ujq). 

Note that in the above when we divide both sides by 
ko the expansion in the Fourier domain is of second or- 
der, where the second order contribution comes from to 2 
term. In our further analysis we will ignore this disper- 
sive term by dropping this second order term. Thus, we 
have the following slowly- varying governing equations for 
the system: 



d_ 

dz 



) Id 

. dz v g dt 



1 d +r 

v dt 
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9- 



(19) 



It follows from equations (19) that the magnetic field is 
proportional to the electric field. Therefore, it is sufficient 
to analyze only the equation for electric field. Also, note 
that if the carrier frequency is in the frequency range that 
determines a positive real value for the phase, then k > 
and the imaginary coefficients will all have a sign opposite 
that which that they are for the case when k < 0. 



AMPLIFYING HOST MEDIA: 
SLOWLY- VARYING APPROXIMATION 

As noted in the previous section, the metallic nanores- 
onators are lossy. To compensate for these losses we pro- 
pose that active atoms be embedded in the dielectric host 
medium or into nanostructures. It should be noted that 
here we consider coherent amplification, where the co- 
herent state of the atomic ensemble is preserved during 
transition from an excited to a ground state. The atomic 
inclusions that we consider here will be two-level, that 
is, they have an excited and a ground state. In the case 
when amplifying atoms are included to host dielectric, 
the polarization, V g , is represented as polarizability, P at , 
of a two- level single multiplied by atomic density, n at : 
'Pg = n a tP a t- The Bloch equations, which describe the 
interaction of polarizability and population difference, N, 
of two-level atoms with the electric field is given in the 
slowly- varying envelope approximation as: 



dN 
~dt 

dVat 

dt 



if 
h 



N£ - i(uj 12 - LO )V at , 



(20) 
(21) 



where loq is the carrier frequency, to\2 is the frequency 
corresponding to the transition between the excited and 
ground state, d is the dipole moment. 

We consider the case when the spectral pulse width, 
Aoj, is within the negative index frequency, Aojni, do- 
main and that the pulse width, To ~ (Aw) -1 , is smaller 
than the characteristic time of the atomic polarization 
relaxation, T 2 . In an experimental situation, AXni <~ 
50nra [13, 14], the value of Ti varies in a broad range. 
For Er ions in a glass matrix at room temperature, 
T 2 - 250/sec and T 2 - 3.6^sec for LaF 3 : PR 3 + at 
2° K [15]. We consider the situation when all these re- 
quirements arc fulfilled and coherent amplification takes 
place. 

To fulfill these requirements and select an optimal 
value for the carrier frequency we analyze the dispersion 
relation (16). The negative index property is determined 
by the magnetic response of the metamatcrial. The pa- 
rameter b determines the strength of this response and 
therefore controls the negative index property. This pa- 
rameter depends on the design of the nanostructure. The 
case of a double split ring was considered in [16, 17], and 
from these papers it follows that the value of b can be 
chosen from the interval [0.1 — 0.7]. In real experimental 
situations of doubly resonant metamaterials, the values 
of ojd and ujt are close. In the case of a fishnet structure 
there is no resonance in the electric response and ujd can 
be chosen to be zero ujd = [16]. In our study we chose 
b = 0.35 and all the frequencies are normalized relative 
to lot- Thus we chose ub = 1-1, lot = 1, w p = 4 and 
the dissipative constants were chosen as 7 = 0.003 and 
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where 



FIG. 1: The dispersion relation 



a = 0.03. This choice of parameters approximately cor- 
responds to the case of silver. The value of 7 is taken 
from [17] and the value of a was extracted from experi- 
mental results presented in [13, 14]. The functions ko(u>) 
and r(w) are presented in Figure 1. The function ko(u>) 
is negative on the interval [1.1 — 1.3]. The function T(uj) 
is positive, which corresponds to the presence of losses, 
and has a flat local minima around to ~ 1.2. The neigh- 
borhood of this flat local minima of T(u>) and interval of 
negative values for ko(u>) overlap. If the carrier frequency, 
luq, is chosen near this local minima then the effect of dis- 
persion of dissipation can be neglected. Since fco(wo) < 0, 
the phase velocity is negative, while on the other hand 
the group velocity, v g , is positive (dko/doj)^ > for this 
interval of frequencies. Therefore, the phase and group 
velocities are oppositely directed for uj selected as de- 
scribed above. If the resonant frequency of the active 
atoms is close to such an uj , the material possesses the 
negative index property and the dispersion of losses can 
be neglected. 



Analytic Computation of the Front 

The balance between coherent amplification and dis- 
sipation leads to the formation of a stationary optical 
pulse. The shape and velocity of such pulse can be found 
using technique suggested in [8, 9]. In the following di- 
mensionless variables: £ = E^e, V g = Pop, N = Nqu, 
C = z/Coj V = t/r/o we can obtain the following governing 
equations for the system 



e-C + e-ri + 7e 

Pr, + i^P 



(1 + iS)p 
-^(e*p + ep*) 
ne, 



(22) 



7 = r oCo, , S 



Co = VOVg, 



Po = No, Eq 



h P 
4??o 



chRc( y/ ii(jj )e(ujo)) 
(27tw Tq v g Re(iJ,(u>o))) ' 



By examining the behavior of the system near the ex- 
cited state, n = 1, we will be able to study the linear 
behavior of the front. The system of equations near the 
excited state has the following approximation in the co- 
moving system of coordinates r = 77 — (/v, z = (: 



(3e r + e z + -ye = (1 + iS)p 
p T + iXp = e. 



(23) 



where (3 — (1 — v)/v. Taking the Fourier transform of 
this system gives us the dispersion relation: 



k(u>) = —/3uj — 27 + 



1 + iS 
uj + A 



(24) 



Thus, near the excited state an amplitude of the electric 
field can be approximated as 



/oo 
e (uj) exp (iujT — i<d(uj)z)duj, 
-00 



(25) 



where Q(u>) — —[3u — i-f + (1 + iS)/(u) + A) is the phase 
of the integral. We are interested in the behavior of the 
traveling wave solution for large distances so as to deter- 
mine solitary wave solutions. The method for determin- 
ing this type of a solution for integrals of this form is the 
stationary phase technique. The phase, ®(ui), has a fixed 
point at: -(1 + i6)/(3 = uj 2 . That is, 



-X±i 




(26) 



For these values of ui we see that the phase has the fol- 
lowing forms: 



Q(oj±) = A/3 - n T 2Wl + i5^/]3 



(27) 



Thus, we see that when the positive square root is cho- 
sen the solution will grow exponentially for large z, which 
cannot occur, and thus the contribution from this term 
must be zero. For the second stationary point there is 
an exponential decay, which implies that we will have 
a strong dependence on z. In our traveling wave solu- 
tion, the shape of the pulse for large z should be de- 
termined by r, and therefore this strong dependence of z 
must be eliminated. Thus, we need the following to hold: 
2Im[iVl + «<5]\//3-7 = 0, that is, 



He 



r 



V fr = 



4:Im{iVTTi8} 2 
4Im[i\/I + i6] 2 
4Im[iVTTi5} 2 + 7 2 
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Using this expression for f3 cr and substituting in lj + 
into the phase of the integral in expression (25) we see 
that near the excited state the electric field can be ap- 
proximated as follows: 

e(z, t) ss eo(w+) exp (iu + r — iRe[Q(uj + )]z)I(z), (28) 

where I(z) is the integral of higher order terms from the 
expansion of the phase around that stationary point w+ . 
Note that Re[6(w+)] = \f3 cr + 2Rc[Wl + idy/j3^). Let 
us use the following notation to simplify this expression: 
VI + iS = S' + iS", which allows us to simplify this term: 

Re[0(a; + )] = ^ ■ (29) 

Also note that if 5" = 0, i.e. if the coefficient in front 
of the polarization in the first equation in the system of 
governing equations is not complex, and if we ignore de- 
tuning, A = 0, then the phase will be zero. Also, note 
that the form of the front depends on e (uj + ), which is 
determined by the initial condition for the electric field. 
In addition, we sec that the shape of the pulse is de- 
termined by the variable r, and the imaginary part of 
w+, that is the shape of the pulse is proportional to 
exp(-Im[o; + ]r) = cx P (-(2( ( 5') 2 /7)t). 

We also see that because lu+ is complex this solution 
will oscillate in r with period determined by Re[0(w+)]. 
Thus, as we have noted before, when 8 = the phase is 
zero, and no oscillations exist. When numerical compu- 
tations are performed this can be observed. 

We have obtained two very important pieces of infor- 
mation from the linear analysis. Firstly, we have deter- 
mined the velocity, v/ r , of the front, and for this reason 
we are now able to solve the ordinary differential equa- 
tion determined by moving into a co-moving frame of 
coordinates. Secondly, we now know that the front os- 
cillates as it propagates, so that before we move into the 
co-moving frame we must express the electric field and 
polarization envelopes as: e(r, z) = a(r, z) exp(iKz) and 
p(r, z) — r(r, z) exp(iKz). Plugging these into Eqs. (22) 
gives us the following system in the co-moving frame 
t = r — z/v 

(-p cr )a t + (iK + 'y)a=(l + i8)r (30) 

n t = ~(ra* +r*a) (31) 

r t + iXr = na (32) 

Near n = 1 we can observe that: 

(~(3 cr )a t + (iK + 1 )a=(l + i5)r (33) 
r t + iXr = a. (34) 

With the velocity determined by our previous analysis 
and recall that K = -(A 7 2 - ^8"8')/A(8'f the charac- 
teristic equation for the above system is degenerate. To 



see this, let us look at the simpler case when A = 0, that 
is, when the carrier frequency is equal to the transition 
frequency of the atomic inclusions. For this case we have 
the following ODE 

- ^a tt + (ij8"8' + 7 (<5') 2 + 7 (<n 2 H - 

(S') 2 ((S') 2 + 2iS'8" - (8") 2 ))a = 0, (35) 

which is degenerate with solutions of the form: 

a(r) = a exp(fcir) + C2rexp(fcir), (36) 

h - 2 « y > 2 + <W "> (37) 
7 

Note that to solve the nonlinear ODE we need both a(0) 
and r(0). If we can determine a'(0) from the initial con- 
dition given for a, then we know r(0). Unfortunately, 
the above linear analysis does not allow us to do this be- 
cause the characteristic equation is degenerate and our 
solution near the excited state does not allow us to sim- 
ply express a'(0) in terms of the o(0). However, if we look 
at the linearization around n = — 1, the ground state, we 
see that: 

(-Pcr)a' + (iK + 7 )o = r (38) 
r' + iXr = -a, (39) 

and using the velocity computed in the previous analy- 
sis our characteristic equation tells us that a(r) can be 
expressed as: 

a(r) = C\ cxp(« + r) + Ci exp(K_r), (40) 

where 

2(1±V2)(£') 2 . 2W(1±V2) 
K± = j + ^ W 

Since we are analyzing the behavior of the front near 
the ground state we know that the shape of the pulse 
must decay exponentially here, and therefore C\ = 0. 
Thus, 

a(r) = a(0) cxp(k_t), a'(r) = k_o(0) exp(«_r) (42) 

Thus, given a choice for a(0) we see that: 

r(0) = ((-&r)«- + {iK + 7))a(0) (43) 

and we are able to numerically solve the nonlinear ordi- 
nary differential equation. These solutions can be seen 
in Figs. (2)-(3) 

To match the asymptotic behavior let us write the form 
of the solution near the excited state in the form: 

a(r,z) = e (r - t„„) exp (fc_ (r - T un ))cxpiKz (44) 

This asymptotic form must match the condition obtained 
using the stationary phase technique. Since the position 
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FIG. 2: This figure represents the numerical integration of 
the nonlinear ordinary differential equation. This corresponds 
to the electric field pulse generated by the transition of our 
initially excited system. 

i.o - 

0.5 - / 



-4 -2 I 2 4 



-1.0 - 

FIG. 3: This figure represents the relative atomic population 
as computed by numerical integration of the nonlinear ordi- 
nary differential equation. Note that as the electric field pulse 
is generated the system goes from the excited to the ground 
state. 

of the solution described by the stationary phase tech- 
nique is described through z while the position of the 
front solution given by the ODE solution is determined 
through the shift variable To, by comparing the asymp- 
totic of these two equations we will be able to express To 
as a function of z. 

Thus, if we add a shift term, To to the ODE solution 
we can compare asymptotic: 

eo(w + ) exp (fc_r + iKz) I(z) = 
eo(r - r - T un ) exp (&-(t - t - T un )) exp (iKz) (45) 

Thus, 

e (w+)exp(fc_T)/(z) = 
eo{T - t - Tun) exp (fc-(r - r - T un )). (46) 

The position of the front is given by To while the width of 
the pulse is given by t— T un . Since the transition from the 
stable to unstable states takes place on the order O(l) 
this implies that t — T un is of order 0(1). On the other 
hand, To the position of the pulse is dependent on the size 
of the incident pulse, eo, that is tq ~ 0(l/eo). Since the 



initialization is small this means that To is large. Thus, 
t >> t - T un and 

e (w+)I(z) = e T exp (fc_T ). (47) 

From this we see that 

ro = (ln(eo(w+)/(«)) - ln(Toeo)) , (48) 

and To, the shift variable, can be determined iteratively. 

Numerical Computation of Nonlinear PDE 

Let us return to the dimensionless system of partial 
differential equations (22). To solve this system numeri- 
cally we use a Strang splitting method, where wc have a 
step that solves the ordinary differential equations, and 
we have a separate step that solves the partial differential 
dimensionless Maxwell's equation. The ordinary differen- 
tial equation step is solved using midpoint method, and 
the partial differential equation step is solved using the 
MacCormack method. At the boundary we consider a 
small Gaussian pulse e(t, 0) = 0.01 exp (— i 2 ) and we also 
assume that the system is initially in the excited state, 
n(0,z) = 1. The result of the numerical computations 
can be seen in Fig. (4). 

As was noted in the section on linear analysis and 
in the previous section, the pulse should oscillate as in 
moves through the medium. In the numerical computa- 
tions this does occur, and the period of these oscillations 
if given by the phase of the approximation of the front 
given in the expression (28). A comparison of the velocity 
computed analytically and the velocity of the pulse gen- 
erated by the system numerically can be seen in Fig. (5). 



THE CASE OF AMPLIFYING 
NANOSTRUCTURES 

Another approach toward loss compensation in mate- 
rials with embedded nanostructures is to add the amplifi- 
cation capability directly into the nanostructures. There 
are two main motivations in favor of this approach. First, 
the small size of the nanostructures will induce highly in- 
tense local electric fields that result from the interaction 
of the nanostructures with an external electromagnetic 
field. This strong local field will enable a strong response 
from the gain material if it is inserted directly into a 
nanostructure. Secondly, since the nanostructures are 
the source of loss, the losses in the system occur on the 
scale of a nanostructure. In the case that was considered 
in the previous subsection, the loss was compensated by 
a gain material distributed over the entire host medium. 
Thus, a strong loss occurred on the scale of a nanos- 
tructure, while the compensatation acted on a mean-field 
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FIG. 4: This figure showed a plot of the intensity of the pulse. 
Analytic and Numerical velocities when complex parts are zero 




0.2 0.4 0.6 



FIG. 5: This figure shows the result of comparing the velocity 
of the front generated by the numerical computation of the 
nonlinear PDE and the velocity of the front that was obtained 
analytically. 



scale. The scales of loss and loss compensation are not in 
agreement when performed in this way. A more appropri- 
ate way to compensate for loss would not only balance 
the gain and loss of energy, but would also match the 
scales of these two procedures. We propose that this can 
be achieved by embedding a gain material directly into 
the nanostructures, and in this subsection we focus on 
modeling such a method of loss compensation. 

In our analysis we consider the simplest case where 
nanostructure in a leading order can be viewed as an 
LC circuit. We consider nanostructures that produce 
spatially homogeneous effective local fields. This simpli- 
fication allows us to avoid a consideration of the spatial 



dynamics for a local field, and to take into account only 
temporal dynamics. The external electric field generates 
a current in the LC circuit. The current in turn estab- 
lishes a local electric field across a "capacitor". Thus, if 
we place active atoms or quantum dots in the part of the 
structure that is responsible for the electric capacity, we 
will be able to utilize the strong local electric field and 
directly compensate for loss. We first derive an equation 
that describes the oscillatory processes in an LC circuit 
that contains active atoms in the capacitor. In the ideal 
case without losses, the potential difference at the capac- 
itor can be expressed in the form: 



U c = -U L +£ V , 



(49) 



Here Ul is potential difference on the inductance and 
£i n is the electromotive force due to magnetic flux varia- 
tions through the area of the LC circuit. We can express 
Ul via the electric field, E c , and the polarization of the 
medium, P a , inside the capacitor: 

U c = d(E c -4iTP a ). 

Here d is the size of the gap in a LC circuit, E c is the 
electric field induced by charges at the edges of a LC 
circuit, and P a is the polarization of the medium inside 
the capacitor. The electric field in the capacitor is 



Ec = 2ES r = Attct = 4tt ( 

. Sc 



where S c is the area of the metallic plates formed at the 
end of the LC circuit. 

Taking into account the expression for E, we see that 
the expression for Uq has the following form: 



U c = d4n (a - P a ) = And ( -|- - P a 



The potential difference on the inductance and the 
electromotive force can be expressed as: 

- LdI f - 

~ c 2 dt ' m " c dt ' 

which means that the equation describing this iC-circuit 
reads 



■4 .-,-</ ( ^ 



L dl 
~&~dt 



that is, 



dl 



Airdc 2 
L 



P, 



1 d<$> 

c dt 



c d$ 
L~dT 



Recall that the capacity is given by: C — Sc/iAnd), 
which allows us to make the following substitution: 



dJ 
~dt 



CL 



(q - S c P a ) 



c d® 
~L~dt 
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From the expression for Thompson's frequency, Wj = 
(LC/c 2 ) -1 , we obtain system of equations 



dl 



c d$ 



dt 



L dt 



dq = 



If we consider a system with loss we add the potential 
drop due to resistance, Ur = IR, to the system and the 
governing equations will read: 



dl R n , . 



7ra 2 c diJ 

" — E d7 



dq = j 

dt 



(50) 



We know that the polarization effects due to the presence 
of active atomic inclusions has the form: 

P a = n a (p* 12 d 12 + pi 2 d 2 i), 

where n a is the density of the atomic population, p 21 
is an off-diagonal element of the density matrix and di 2 
is the dipole moment. Since these atoms are contained 
within the capacitor, we can express the atomic density 
as: 



n a = 



N, 



S c dJ 



where E is the external electric field, P ns is the polariza- 
tion contribution from the nanostructures, 5 is the losses 
due to plasmonic oscillations, u>d is the frequency of the 
dimensional quantization that comes from the geometry 
of the resonator, and u> p is the plasma frequency. 

Now that we have the material equations for this case 
we need to derive the slowly- varying envelope equations. 
From Eq. (51), we take the Fourier transform to see that: 



M = 



pu 2 



-H 



(53) 



Thus, 



B 



1 + 



4tt/3w 2 



H 



lA ^ a -P a ,(54) 



- UT 



and we see that the effective permeability of the system, 
jj,(u>), is the same as in the case of mean-field amplifi- 
cation. By taking the Fourier trasform of Eq. (52) we 
obtain 



D 



J D 



iSuj 



E 



(55) 



to see that the effective permittivity, e(u), is also the 
same as that for mean-field amplification, which is de- 
termined by the nano-structures. Recalling the Fourier 
transform of Maxwell's equations we have: 



where N a is the population and S c d is the volume of the 
capacitor. Thus, the polarization term that comes from 
atomic inclusions placed inside the capacitor of the LC 
circuit can be expressed as: 



k 2 E 



E = 



" Pns ~t~ 



M 



k 2 H 



L0 2 fl Q (uj) 



H 



4ttuj 2 £o(lu) ~ Amok ~ 
m + P ns 



N a 

Sc'Pa = -f(Pl2 rf 12 + Pl 2 d 2 l). 

a 

Since the size of the "capacitor" is small it follows that 
the polarization contribution is quite large. Recalling 
that the magnetization, M, can be expressed as a func- 
tion of the current in the LC-circuit, M — xl{t). The 
material equation determined by these nanostructures 
with inclusions placed inside the capacitor is therefore: 



d 2 M n 
dt 2 



+ 7 



dM ns 
dt 



d 2 H 



dt 2 +UT°- Br , (51) 



where b = 'Ka 2 c/Ln, a — n/d, and 7 = R/Lk. 

Now, as in the previous case we will be consider- 
ing structures responding as an iC-circuit to the mag- 
netic field and these structures can have a resonant re- 
sponse to external electric field. The polarization effects 
comes from plasmonic oscillations and the geometry of 
the nanostructures: 



where M comes from the nanostructures. Note that the 
polarization effects from the active atoms interact mainly 
with the local electric field inside the capacitor, which is 
much stronger than the external electric field. Thus, we 
have 



2 iv 2 i(uj)fi(uj) 



k 2 - 



E 



H 



i4:TTu> 2 kLO^a 
c(wj, — uj 2 — iivy) 
i(4w) 2 uj 3 e (Lu)uj 2 ,a 



^ (56) 



c 2 (u 2 



3(57) 



Thus, the expansion of the bracketed terms, and the 
eventual resulting slowly-varying envelope approxima- 
tion is similar to the previous case. 



d_ 

dx 



— + — — + r 

dx v g dt 
$ d 



Vg dt 



in 



s = - 

H = - 



iATrojQkoco^a 

i(47r) 2 WQeo(wo)w 2 -o: 
c 2 (uj, — oJq — iivoj) 



d 2 P n; 

dt 2 



dt 



- LO D P ns 



4^' 



In addition to the above material equations we must re- 

(52) 

call that the time derivative of the local electric field is 
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directly proportional to the magnetization that results 
from the nanostructures, that is, 



dE, 



loc 



dt 



= vM n 



47T 

S r .K 



(60) 



Thus, we have that the Fourier transform of the local 
electric field can be written as the direct sum of the 
Fourier transforms of the external electric field and po- 
larization contribution from the active atomic inclusions, 
that is, 



E, 



loc 



-H + 



Pa (61) 



where as before H is the external magnetic field. 

We now add to these slowly- varying equations those for 
the polarization and relative atomic population as deter- 
mined in the previous case. The difference we must keep 
in mind is that these amplifying structures arc embedded 
into the nano-structures, and therefore the electric field 
in these equations is the local electric field rather than 
the external. The slowly-varying governing equations for 
the amplifying material: 

-^T = -rN£ loc - l{u}\2 - LU )Va 

dt h 
dN i 

~~Qj~ = ~ J^i^locPa + £loc'Pa)i 

where as before a; 12 corresponds to the transition fre- 
quency and ujq is carrier frequency. Thus, W12 — u)o is 
the detuning. To simplify this expression let us replace 

&loc * i&loc- 

Thus the slowly-varying envelope approximation for 
the governing equations are: 



— + — — + r 

dx v g dt 



dV a 



i(47r) 2 a>g6 (^o)t4a 

n = - 9 , 2 2 — : — \ ' a y b2j > 

C^(Ll>t — Ll!q — iujoj) 

'-N(rjH + XV a ) - i{u; 12 - uj )T a (63) 



dt 
that is 



dt 

:1 V = ~((rfH + XP a )*V a + ( V H + XV a )V*), (64) 



d_ J_d_ 

dx v g dt 



7i = aV a (65) 



<¥P 2 

—^ = -( 11 NH + \NP a )-i(u l2 -u a )P a (66) 
' V =-U v *H*P a + r 1 HP* a + (X* +X)\P a \ 2 ), (67) 



dt 



where a is i times the complex coefficient in front of the 
polarization from the previous equations. This is done to 
make notation simpler. 

To determine the governing equations in dimensionless 
variables let us first multiply the first equation in this 



system of equations by 77 and let us relabel -qH = H for 
convenience. Thus we will have: 



— + ——+T 

dx v„ dt 



H = aV a (68) 



OVa 



■(NH + XNVa) - »(W12 - u )V a (69) 



dt h y 

,,V = -±. iH * Va+HK + { \* +X )\P a \ 2 ), (70) 



dt 



where a = r\a. Now, let us write a = a' + ia" and in- 
troduce the following dimensionless variables: H = Hoh, 
P a = P p, N — N n , x — z/zo, t = t/r a . Thus in 
dimensionless variables the governing equations read 



h z + h T + 7/1 = (1 + ia)p (71) 
p T = nh + Xnp + iAup (72) 

-l(h*p+ hp* + (X* + X)\p\ 2 ) (73) 



n T = 



where z = T v g , P = N , 7 = Tz , X = X/(T v g a'), 
a = a" /a', Alo = tq(u — w 12 ). This is a new sys- 
tem of equations describing electromagnetic field inter- 
action with nanostructures and the interaction of nanos- 
tructures with gain material via the induced local fields. 
Note that the usual conservation law holds in this case: 
|p| 2 +n 2 = l. 



Front velocity computation 

To compute velocity of a steady wave solution we will 
use technique described in Near the excited state n w 1. 
This reduces the system of equations to 



h x + h t + = (1 + ia)p 
p T = h + Xp — iAuip 

In a co-moving frame r = t — x/v, z = x 

h z — (3h T + 7/1 = (1 + ia)p 
p T = h + Xp — iAujp 



(74) 
(75) 



(76) 
(77) 



where (3 = (1 — v)/v. 

The dispersion relation for this equation therefore 
reads: 



k(u) — —(3uj + i-f + 



1 + ia 



oj — Au) + iX 



(78) 



We know from the previous case that to obtain the ve- 
locity of the front we must find the zeros of the first 
derivative of the dispersion, which in this case are: 



uo± = Auj + iX±i 




(79) 
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Using these to evaluate the imaginary part of the phase 
as before we obtain: 

Im[fc(w±)] = -/3X' + 7 =f 2Rc[iVl + i&] \Tp (80) 

where A = A' + iX". Thus, the velocity is given 



-2Re[i>/i + ia] + W^Re^vTTla]) 2 + 4A't~ 
VP = ;ru ( 81 ) 



2A' 

If we expand around A' = we see that up to second 
order 



[3; 



[4; 



[5; 



7 



7 2 A' 



7 3 (A') 2 



+ ■ 



+ ia 



From this we see that when A' — the velocity is the 
same as the most basic case if a = as well, that is 
(3 = 7 2 /4. 



CONCLUSION 

In this paper we theoretically studied the propagation 
of a steady state ultrashort pulse in a dissipative medium. 
We considered the case when the medium consists of lossy 
metallic nanostructures embedded to a gain material and 
the case when the gain material is embedded directly to a 
nanostructures. We derived a system of governing equa- 
tions and found analytic expressions for the velocity and 
spatial period of an optical pulse coupled with the polar- 
ization wave. We also proposed a simple technique that 
describes the shape of this pulse using a simple system 
of ordinary differential equations. 
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